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[JLh : ABSTRACT 

CN I A new algorithm for implementing the adaptive Monte Carlo method is given. It is 

T-H ■ used to solve the Boltzmann equations that describe the time evolution of a nonequilibrium 

C^^ ' electron-positron pair plasma containing high-energy photons. These are coupled nonlinear 

oo : 

T-H . integro-differential equations. The collision kernels for the photons as well as pairs are eval- 

(N ■ 

O ' uated for Compton scattering, pair annihilation and creation, bremsstrahlung, and Coulomb 

i> : 

Q\ ■ collisions. They are given as multidimensional integrals which are valid for all energies. For 

'"5 ' an homogeneous and isotropic plasma with no particle escape, the equilibrium solution is 

Q . expressed analytically in terms of the initial conditions. For two specific cases, for which the 

$—1 ' 

^ [ photon and the pair spectra are initially constant or have a power law distribution within 
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the given limits, the time evolution of the plasma is analyzed using the new method. The 
final spectra are found to be in a good agreement with the analytical solutions. The new al- 
gorithm is faster than the Monte Carlo scheme based on uniform sampling and more flexible 
than the numerical methods used in the past, which do not involve Monte Carlo sampling. 
It is also found to be very stable. Some astrophysical applications of this technique are 
discussed. 
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1 Introduction 

Nonthermal emission of liigh-energy radiation from a variety of compact astrophysical 
objects e.g., 7-ray-burst sources (Meszaros & Rees 1993a,b), pulsars (Chen & Ruderman 
1993), active galactic nuclei (AGN; Lightman & Zdziarski 1987; Svensson 1994; and Padovani 
1996), and jets in the AGN (Sikora 1994) seem to indicate the presence of a relativistic 
electron-positron pair plasma in the dense radiation fields of those sources. Such plasmas 
may exist also in the accretion disc coronas of the Galactic X-ray binaries (Sunyaev et al. 
1992), the ergo-spheres of Kerr black holes (Piran & Shaham 1977), and the black- hole 
accretion discs (Tanaka & Kusunose 1985; Gunnlaugur & Svensson 1992). It is conceivable 
that the pair plasma in some of these sources is in thermodynamic equilibrium with itself and 
probably in equilibrium with the radiation. However, it is more plausible that many of them 
may consist of nonequilibrium pair plasmas (Coppi & Blandford 1990 - CB90 henceforth; 
Zdziarski 1988 and 1989). Many of the previous papers on this topic have concentrated on the 
the properties of a relativistic pair plasma in thermal equilibrium (e.g., Bisnovatyi-Kogan, 
Zel'dovich, & Sunyaev 1971; Lightman & Band 1981; Lightman 1981 and 1982; Svensson 
1982b - henceforth S82b; and Zdziarski 1985). Examples of the time evolution of a thermal 
pair plasma, taking into account the finite-medium radiative transfer effects can be found 
in Guilbert & Stepney (1985) and Kusunose (1987). There are not many papers that deal 
with the evolution of a nonequilibrium pair plasma in detail; some examples can be found 
in Lightman & Zdziarski (1987), Svensson (1987), Zdziarski, Coppi, & Lamb (1990), and 
Coppi (1992 - C92 henceforth). 

These investigations are generally based on the Monte Carlo (MC) methods or on 
solving the Boltzmann equations (kinetic theory approach). In a simple MC method based 
on uniform sampling (Pozdnyakov, Sobol' & Sunyaev 1977), individual particles are followed 
as they undergo interactions in the source. In this method, it is usually easy to take into 
account the spatial inhomogeneities and radiative transfer effects well. But it typically suffers 
from relatively poor photon statistics at higher energies and does not lend itself to time- 
evolution calculations involving broad-band spectra. For examples of such MC simulations, 
see Novikov & Stern (1986). In the kinetic theory approach, the system is represented by 
the photon and particle distribution functions which are discretized in energy as well as the 



spatial coordinates and the time evolution is determined by solving the Boltzmann equations 
numerically. In general, it is very difficult to solve the resulting integro-differential equations. 
Moreover, they are usually "stiff" (i.e., there are very different time-scales in the problem). 
The principal advantage of this approach is that it gives good photon statistics at higher 
energies. Some examples of the kinetic theory approach can be found in C92, Ghisellini 
(1987), Svensson (1987), and Fabian et al. (1986). 

There have been some attempts to improve the photon statistics in the conventional 
MC schemes which go by the name phase-space density (PSD) array representation. In this 
approach, the system is represented by the discretized distribution functions (as in the kinetic 
theory approach) but the particle or photon transitions between the energy bins is simulated 
using the MC method and the interaction between the spatial cells is modeled with the aid 
of the escape probabilities. So far this approach has been used to model only homogeneous 
and spherically symmetric systems (e.g.. Stern 1985). Another recent variant of the MC 
method is based on the large-particle (LP) representation, which is described in detail by 
Stern et al. (1995). In this scheme, the system is represented by an array of "large particles" , 
each of which corresponds to a group of real particles sharing the same physical parameters 
(i.e., particle type, position, momentum, and energy). It is more flexible than the PSD 
approach in the sense that each LP is tagged with a statistical weight, which is proportional 
to the number of real particles represented by that LP. For example, this weight can be 
assigned based on the total energy carried by each LP. In many nonequilibrium systems of 
interest in astrophysics, the number of particles in the low-energy range is typically several 
orders of magnitude larger than that of the particles in the high-energy range. Therefore, 
the efficiency of the method may be improved by assigning lower statistical weight to the 
low-energy LPs. Intuitively this approach makes sense but there is no general proof for 
its validity or effectiveness (except for the numerical experiments presented by Stern et al. 
1995). Besides, the statistical weights are rather ad hoc. 

From the preceding discussion, it is clear that the main problem in the analysis of 
nonequilibrium pair plasmas is the computational difficulty. The principal aim of this paper is 
to present a new method for solving the kinetic equations based on an adaptive MC sampling 
scheme. It is faster than the conventional MC method (based on uniform sampling) and is 



more flexible (and in some cases, faster) than the numerical methods previously used. Our 
method resembles the LP method described above, in the usage of the statistical weights, 
but it is much more rigorous. Moreover, it can accommodate anisotropic distributions with 
greater ease. 

In a relativistic plasma containing arbitrary densities of pairs and the high energy pho- 
tons, the collision cross sections for various microscopic processes depend on the energy. One 
cannot use, for example, the simple Thomson cross section as one can do in the nonrela- 
tivistic case. In addition there is a creation and annihilation of the pairs and photons that 
alter the densities. Therefore we have to follow the time evolution of the number density 
as well as the spectrum of each species. Besides, the problem is inherently nonlinear due 
to the form of the collision kernels in the Boltzmann equations. It is possible to write all 
the collision kernels as multidimensional integrals. This reduces the problem of solving the 
coupled Boltzmann equations for the photons and the pairs into a purely computational task 
of evaluating many of these integrals, after each time step, quickly and efficiently. This way 
of formulating the problem of kinetic theory is more ffexible in accommodating any kind of 
distribution functions. We have developed a new algorithm, based on Monte Carlo sampling, 
for computing such integrals. The functional form of the integrands is not assumed a priori. 
Also, no constraint is placed on the shape of the integration region. Usually such integrals are 
evaluated either numerically (by using an equally spaced discrete grid) or through a Monte 
Carlo sampling technique. In order to make the former method more efficient, we have to 
choose the shape of the discrete mesh depending on the form of the integrand. This takes 
away the ffexibility from the method (i.e., the algorithm will depend on the form of the in- 
tegrand). The latter method, based on uniform sampling throughout the integration region, 
is widely used in astrophysics. It is possible to speed up the computation in this method, 
by sampling selectively i.e., sampling more frequently in those domains where the integrand 
is larger. This scheme is known as the importance sampling method or the adaptive Monte 
Carlo method. There is an algorithm, originally due to Lepage (1978), which implements 
this. However it is not well suited for the type of integrals that arise in the present context. 
We have developed a new algorithm to implement the adaptive Monte Carlo method which 
is very efficient (see below). 



In the next section we define various quantities, explain the basic pair plasma model 
we use, and write down the general kinetic equations. In sections 3 and 4 we give the 
integral expressions for various collision kernels, that are valid for all energies. These collision 
integrals are cast in a form that is well suited for the Monte Carlo integration. In section 
5 we describe how we integrate the Boltzmann equations numerically. There we explain 
the adaptive Monte Carlo algorithm we use. In section 6 we express the final equilibrium 
state of an homogeneous and isotropic plasma (with no escape of particles or photons) 
analytically in terms of the initial conditions. Then we apply our time-evolution code to two 
specific examples of nonequilibrium configurations and compare the final results with the 
corresponding analytical solutions. These examples serve as a test for the overall formalism 
presented in this paper. Finally, in section 7, we summarize this work and discuss some 
astrophysical applications. 

2 Model, definitions, and the notation 

We consider a neutral, stationary, and unmagnetized pair plasma which is nonthermal 
(i.e., not in equilibrium). We assume that the plasma is homogeneous and isotropic. If 
the plasma is in a moving source we must interpret all the physical quantities given below 
as the comoving- frame quantities. The number densities (i.e., the number of particles per 
unit volume) of the electrons, positrons, photons, and protons are given by: n_,n^,n^, and 
Up, respectively (n_ = n_|_ + rip). Throughout this paper we express the momentum and 
energy in units of mc and mc^, respectively. Here m is the electron rest mass and c is the 
speed of light in free space. Therefore the momenta and the energies of the particles, as well 
as the photons, are represented by dimensionless numbers everywhere. For the models we 
consider here the protons can be assumed to be at rest. We assume that the state of the 
plasma is completely described by the Lorentz invariant distribution functions f±{x,p) and 
f^{x,p), for positrons, electrons, and photons, respectively. Here x,p represent the position 
and the momentum four-vectors, respectively and x, p represent the corresponding three- 
vectors. Our choice of the metric is such that p^ = 1 for electrons. In the case of photons 
we have p = £:(!, k), where e is the photon energy and k is a unit vector in the direction of 
its three-momentum. Similarly, p = 7(1, /3) for the pairs. Here 7 is the Lorentz factor and 
f3 is the velocity in units of c. We denote the magnitude of /3 by (3. The number density of 



the particles of type i with a momentum p is given by fi c/^p. We define the total densities 
of various species to be rii = J fi{p) d^p, where the integration extends over all values of 
the momenta. Because of the isotropy, we have d^p = Ane'^de in the case of photons and 
d^p = An (3 7^ d'y for the pairs. 

Since we assume that the plasma is homogeneous and isotropic, various distribution 
functions depend only on time and the energy (or the magnitude of the momentum). We 
define the spectral functions for photons, positrons, and electrons to be 

F,{e) = f,{e) and ^±(7) = ^^/±(7), (2.1) 

respectively. The time dependence of these functions is not shown explicitly. The spectral 
functions are normalized so that 

poo /"OO 

/ deF^{e) = 1 and / d-^F±{-f) = 1. (2.2) 

«/ (J «/ 1 

We see that the number of photons of energy e per unit volume and unit energy is given 
by n^F^{e). We will assume that the electrons and the positrons have the same spectral 
functions i.e., -^-(7) = -^+(7) for all values of 7, which we denote by -Fe(7)- 

The equilibrium spectral functions, which are independent of time, are given by 

1 e^ 

^'^'^ ^ 203)0^ exp(£/e) - 1 ^^-^^ 



and 



^=W = ex;(I7e)''-'''-p(-'/e)- (2-4) 



Equation ( |2.3| ) comes from the Planck function for the photons, where ( is the Riemann zeta 
function and C(3) = 1.202. In that equation we have used the equilibrium density of photons 

n^ = 167rC(3) l^—Qj , (2.5) 



where h is the Planck's constant. Equation ( p.4| ) is the relativistic Maxwell-Boltzmann 
distribution for electrons and K2 is the second order modified Bessel function of the second 
kind. In all these equations G = ksT /mc^ is the dimensionless temperature of the plasma, 
where T is the temperature and ks is the Boltzmann constant. 



To study the time evolution of this system we should proceed from the relativistic 
Boltzmann equations for the pairs and photons. In the latter case it is the same as the 
radiative transfer equation. The Boltzmann equation (see e.g., de Groot, van Leeuwen, & 
van Weert 1980) for the particles of type i, described by /j, which takes into account the 
collisions with the particles of type j, described by fj, is given by 

p^d,U{x,p) = Y^f^dn' mx,p')f,ix,q') - Mx,p)f,ix,q)] Fa,,. (2.6) 

Here 9^ is the partial derivative with respect to x^^ and the summation over /i is implied. 
The summation for j extends over all relevant processes. Here g° is the energy component 
of the 4-vector q. Using the initial and the final momenta to designate the particles, the 
collision processes can be represented as p + g ^^ p' + q'. The solid angle around one of 
the outgoing particles is dQ'. Finally, (Tjj is the cross section for the process and F is the 
invariant flux factor. It is necessary to remark that in the present form, the above equation 
cannot account for the quantum mechanical Bose enhancement and Fermi blocking effects, 
respectively for the photons and pairs. In order to do so, we need to take into account the 
particle occupation numbers in the phase space. For photons, this is given by 

ftW = i(-)'/,(^)=f-)'^. (2.7) 

2 \mcj \mcj Sne^ 

which in the equilibrium case reduces to l/[exp(£/0) — 1], as expected. If we are considering a 

process in which two particles of momenta p and q produce a photon of momentum p', then we 

should make the replacement fi{p)fj{q) —* /i(p) /?(?)[! + 9-y{p')] in the Boltzmann equation. 

These effects play a significant role only when (^^ ~ 1 or n^e^'^F^{e) ~ 1.76 x 10^° cm^^. 

For the densities and the energies of interest here, these quantum mechanical effects can be 

neglected. An analogous remark applies to the case of the pairs. Such induced effects in 

a relativistic thermal plasma at high temperatures and densities have been considered by 

many authors in the past (e.g., Ramaty, McKinley, and Jones 1982). 

The Boltzmann equations reduce to simple rate equations in the comoving frame as 
a result of the homogeneity and isotropy of the plasma. We denote the comoving time 
coordinate by t. The rate equations are given by 

^f^ =j:h- h X^], , (2.8) 



where i stands for either photons or electrons and q labels the binary collision process (Comp- 
ton scattering, pair processes, bremsstrahlung, or Coulomb collisions). The summation runs 
over all those processes that involve a particle of type i among the products of the collision. 
Here rji is the emission coefficient for the production of a particle of type i with momentum p 
(or scattering of such a particle into that final momentum state) and Xi is the corresponding 
absorption coefficient. Notice that /j, rji, andxi depend only on the energy of the particles 
and time. In order to obtain the collision kernels, r/j and Xi? we require the binary reaction 
rates in a relativistic plasma (e.g., de Groot, van Leeuwen, & van Weert 1980; Baring 1987a). 
Using the appropriate reaction rates we can write 



f^ l + dim Ju dP 



(2.9) 



where dim = 1 for identical colliding particles (i.e., / = m) and is zero otherwise. The 
summation in this equation is over those incident states (labeled by / and m) which result 
in a final state labeled by i. Furthermore, daim/dP is the differential cross section for the 
process whereas dP is a shorthand for d^p which is defined above. The four-momenta of 
the colliding particles are given by pk = {pi, Pk) (for k = l,m) and the four-momentum of 
one of the outgoing particles is p. The product of the phase-space densities of the colliding 
particles dFim is given by 

dFim= n fAPj)d'Pr (2-10) 

j=l,m 

We have d^pi = efdeidfli for the photons and d^pi = (3i 'jf dfli d'ji for the pairs. The kinematic 
factor J-'im for binary collisions (see e.g.. Landau & Lifshitz 1975) is given by 

^Im = {Ul ■ Um)Pvel{Pl,Pm), (2.11) 

where ui = pi/p^ and /^rei is the relative velocity of the colliding particles in units of c. If at 
least one of the colliding particles is a photon we will have /3rei = 1- Otherwise 

PveliPhPm) = ^ ^ ■ (2.12) 

^ ~ h>r Pm 



The integration in equation (|2.9|) is over a region U of the phase space of the colliding 
particles, which is specified by the energy-momentum conservation. It depends on the energy 
p" of the final state. Now we specialize to the case of a process for which the reacting particles 
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are labeled by / = 1 and m = 2. By using equation (|2.1|) we can express dFi2 in terms of 
the spectral functions and the densities. This gives the following final expression for the 
emission coefficient (i.e., the production rate) for electrons or photons: 

Now we define the total reaction rate between two particles of energies ei and 62 to be 

R(gl,g2) = „., ,W / C^/iJ'iaCTtotal, (2.14) 

2(1 + Ou) J-i 
where /i is the cosine of the angle between the momenta of the colliding particles and o"totai 
is the total cross section for the process considered (integrated over the entire phase space 
of the emitted particle). Clearly J^i2 as well as crtotai depend only on ei, 62, and /i. Now it 
is possible to express the emission coefficient in terms of the total reaction rate as 

V{e) = J ll[de,F,i6,)]R{e,,e2)P{e,,e2;e), (2.15) 

where the integration is over all values of ei and 62 without any restriction (in contrast 
with eg. [p.l3|| ). In the above equation, P is the probability, integrated over all incident and 
emergent angles of the particles, for emitting a particle of energy e, from a collision between 
the particles of energies ei and 62- It is normalized so that J deP{ei,e2',e) = 1, where the 
integration is over all values of e. Equation ( |2.15| ) has been used by several previous authors 
(e.g., CB90). 



We can obtain the absorption coefficient from equation ( p^.9| ) with only minor changes. 
For the absorption of the particles of type i with a momentum pi we find that 

Here atotai is the total scattering cross section for the process. The summation extends over 
all relevant processes. For a binary process, involving the particles of type i and type j, the 
absorption coefficient can be written in terms of the spectral functions as follows: 



where the integration region U is determined by the energy-momentum conservation. This 
way of writing the emission and absorption coefficients is very convenient for Monte Carlo 
evaluation we describe below. 



We remark that in equations (|2.13| — |2TT7| ) we have used e in a generic way and it has 



to be replaced by 7 whenever it refers to the pairs. Physically, Aire'^'qle) is the rate at which 
photons of energy e are emitted per unit volume and unit energy due to the process under 
consideration; similarly, An (3 'y'^ r]{'~f) gives the corresponding electron emission rate (recall 
that we express energy in units of mc^). Electron and photon absorption rates are obtained in 
a similar way. If the size of the system is /, the optical depth r and the absorption coefficient 
are related by r = Ix/c. Equations (|2.13|) and ( p.l7|) constitute the point of departure for 



the following two sections where we obtain the emission and the absorption coefficients for 
the photon and the pair kinetic equations. We remark here that in the case of Compton 
scattering of the photons as well as the pairs, the collision integrals only give the rate at 
which the spectrum changes at a given energy and do not imply any change in the total 
numbers of the particles. 

3 Collision integrals for photons 

The preceding discussion has been very general. We now obtain the integral expressions 
for the photon emission coefficients due to Compton scattering, two-photon pair annihilation, 
and bremsstrahlung and the absorption coefficients due to Compton scattering and the pair 
creation. In this paper we do not consider the double-Compton emission or the three-photon 
emission through pair annihilation . Also we do not consider the effect of photon absorption 
through the inverse-bremsstrahlung (free- free absorption). 
3.1 Compton scattering of photons 

The problem of Comptonization in astrophysics has been analyzed extensively by many 
previous authors (e.g., Blumenthal & Gould 1970; Rybicki & Lightman 1979; and more 
recently by CB90). Here we obtain an integral expression which is valid at all energies of 
the incident electrons and photons. Throughout this paper we call the comoving frame of 
the plasma the C-frame. Let p and pi be the momenta of the incident electron and photon, 
respectively in the C-frame. Let q and gi be the corresponding momenta after the scattering. 
Recall that p^ = 1 and pf = 0. We require the final photon energy to be e. Hence we set 
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gi = e(l, k), where k is the directional unit vector. We write p = 7(1, /3) and pi = £1(1, ki). 

Here 7 is the Lorentz factor of the incident electron, (3 is its three- velocity in units of c, ei is 

the energy of the incident photon, and ki is its directional unit vector. Using the fact that 

{p + Pi — qi)^ = q^ = 1 we obtain the well known relation between the initial and the final 

photon energies viz., e = e{ei) or ei = ei{e), where 

ai7ei ^ aje 

and £1 = — . (3.1) 



aj + bei ai7 — be 

Here a=l — /3-k, ai = l— /3- ki, and 6=1 — cos^, while cos^ = k • ki gives the cosine 
of the photon scattering angle in the C-frame. Let fi = cos 6 and the cosine of the angle 
between /3 and k is defined to be fi'. The angle between the planes formed by the pairs of 
vectors (k, ki) and (k, /3) is defined to be 0. It is easy to show (see the appendix for further 
details) using equation ( p.l3|) that the Compton emissivity for photons is given by 



vie) = "''"^'l" V^^"' / (rf7ci/irf/i'#)F,(7)F,(£~i) (,r^] , (3.2) 

where A = ,^^ — ^ sin^ ^' + 1 and S, = 017/(017 — fee), while 9' is the photon scattering angle in 
the rest frame of the incident electron and rg is the classical radius of an electron. The region 
of integration U is defined by 7niin < 7 < 7max, —1 < Z^, /w' < 1, and < < 27r subject 
to the condition that eimin < £^1 < ^imax- Here 7min and 7max are the limiting electron or 
positron energies in the plasma. Similarly eimin and ^imax are the limiting photon energies. 

Now we obtain the corresponding "absorption" coefficient (as stated before, this is not 
a real absorption; the photons are scattered into a different energy bin). Let p = e(l, k) and 
q = 7(1, /3) be the initial momenta of the photon and the electron, respectively. Various 
symbols have the same meaning as above. The photon energy in the rest frame of the incident 
electron is given hj x = p ■ q = 'ye{l — jS fi), where /i is the cosine of the angle between the 
vectors /3 and k. Now rij = n+ + n_, 6ij = 0, Fj = Fe, dQj = 2Tidji, Tij = (1 — /9/i), ^i = £, 
and Ej = 7. Substituting these expressions into equation ( p.l7| ), we obtain 

= ^i^-+^+) f dftdjF^i^Kl - /3/i)atotai(x), (3.3) 

2 Ju 



where 

^2a;(l + x) 



o-totai(a^) = 27rrg <^ — — 



l + 2a; 



In (1 + 2x) 
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ln(l + 2x) l + 3x ] 

2x (l + 2x)2j ^^ 



is the total cross section for Compton scattering (e.g., Jauch & Rohrlich 1980 - JR80 hence- 
forth). The integration domain U is defined by 7min < 7 < 7max and —1 < /i < 1 without 
any restriction. Here 7min and 7max are the hmiting electron energies, as in the previous case. 

3.2 EMISSION AND ANNIHILATION OF PHOTONS BY THE PAIRS 

The emissivity due to the annihilation of relativistic electron-positron pairs (creating 
two photons) has been analyzed by many authors before (e.g., Zdziarski 1980; Ramaty & 
Meszaros 1981; Yahel & Brinkmann 1981; Svensson 1982a - henceforth S82a). We give here 
the final result using the notation of S82a and refer the reader to that paper for a detailed 
derivation. Let Pi = 7^(1, /3j); i = 1,2 he the momenta of the electron and the positron, 
respectively in the C-frame. Let gi = e{l, k) be the momentum of one of the emitted photons. 
Here c(3^ are the particle velocities and 7j are the corresponding Lorentz factors, e is the 
photon energy, and k is its directional unit vector. The momentum of the C-frame itself is 
denoted by g = (1, 0). We call the center-of-momentum frame of the pair the CM-frame and 
the quantities in this frame appear with a suffix 'cm'. The particle momenta in this frame 

are Picm = 7cm(l, /3cm)) P2cm = 7cm(l, -/3cm)> <?! cm = £^cm(l,kcm), and gem = 7c(l,-/3c)- 

Here 7cm is the Lorentz factor of the electron or positron, ^cm is the photon energy, and k^m 
is its directional unit vector (in the CM-frame). The velocity of the CM-frame as measured 
in the C-frame is c f3^ and 7c is the corresponding Lorentz factor. Various directional cosines 
are defined as follows: /i, x, y, and z are the cosines of the angles between the pairs of vectors 
{/3i, (32), (kcm, /3cm)5 (/^c) /^cm)' and (/3c, kcm), respectively; the angle between the planes 
formed by the pairs of vectors (/3c, kcm) and (/3c, /3cm) is denoted by 0. After analyzing 
the kinematics, we obtain 7cm = ^[1/2 + 7i72(l - /3i/52/^)/2], 7c = (71 +72)/(2 7cm), y = 
(7i-72)/(2/?c/5cm7c7cm), z= (^ - 7c7cm)/(/5c7c7cm) , and X = yz + y/[{l - y"^) (l - z"^)] COS (f). 
Now using equation ( p.l3|) we obtain the pair emissivity 



cn+ri- 



2 



/?cm 7cm ( d(j\ 



vie) = ^-P^ / d^^d<|> UiFeil.) ^7.] -P^^ (^ . (3.5) 

4 7re^ Ju fj^ /5c7c7i72 V^^/cm 



The differential cross section in the CM-frame is given by 



2 



dcr\ r^ 



dn _ 4/3cm7c^, 



cm 



27 



-l + ^r^(C+ + C-)-7r:^ a + C^ 



cm 



(3.6) 



where (± = 1/(1 ± /^cm^^)- The integration domain U in equation (|3.5| ) is given by 7min < 



7i,2 < 7max, — 1 < /^ < 1, and < < 27r, subject to the condition — 1 < 2 < 1, which 
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is equivalent to the condition r_(7i, 72, /i; e) < 7cm (71, 72, /w) < r+(7i,72,/i; e), where T± = 
£:7c(l ± Pc)- Here 7inin and 7max are the hmiting pair energies in the plasma. 

Now we obtain the photon absorption coefficient due to pair creation. Let the initial 
momenta of the photons be p = £:(1, k) and p' = £'(1, k'), with the usual meaning for various 
symbols. If an electron-positron pair is produced, then the CM-frame Lorentz factor of the 
electron is given by 7cm = Vl^^'i^ ^ f^)/'A^ where fi is the cosine of the angle between the 
vectors k and k'. Using equation ( |2.17| ) we find 

CTi f 

X{e) = ^J^ dfide'F^{e'){l - /i)(Ttotai(7cm). (3.7) 

Since (7(77 — > ee) = 2P^^a{ee — > 77), by integrating equation ( ^^ , we find 



0"totall7cmJ — 2 



In , , „ 

Pern V 1 - Am / 7,^ 



cm 



(3. 



leva 

The integration domain U in equation ( p.7| ) is defined by — 1 < /i < 1 and e* < e' < ^max, 
where e* = 2/[e{l — /i)] is the pair creation threshold energy and ^max is the limiting photon 
energy in the plasma. 

3.3 BREMSSTRAHLUNG EMISSIVITY 

The bremsstrahlung emissivity of a pair plasma has been analyzed in several papers 
(e.g., Haug 1975b, 1985c, 1987, 1989 and Dermer 1986). The final expression for the photon 
emissivity can be written as 

%ai.(e) = ^ f dfj^dn n[Fe(7.) ^7.] — \l (< + nl) ^ + n+n.-^l , (3.9) 

where a is the fine structure constant. The ffist term inside the braces represents the sum 

of the electron-electron and the positron-positron contributions and the second term gives 

the electron-positron contribution. The expressions for p, JF12, Cj, and Aj, along with the 

definitions of the integration variables fi and Q are given in the appendix. The emissivity 

due to pair-proton bremsstrahlung can be written as 

cnp(n+-Fn_) nmax ( da\ 

ri,.oUe) = -^^:^2 1 dlF^W[-^^) , (3-10) 

\ / proton 

where ((icr/rfe) proton is the cross section for this process (see e.g., JR80). Here the protons 
are assumed to be at rest. 

4 Collision integrals for pairs 
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Figure 1: The emissivity due to (a) electron-electron and (b) electron-positron 
bremsstrahlung from a thermal plasma for three different temperatures. The dashed lines 
represent the emissivity due to pair-proton bremsstrahlung (given here for comparison) . The 
energy of the emitted photon is e and S = Aire^rile) / {cnin20'Th) , where ni,2 are the appro- 
priate densities. These results agree with Haug (1985c) and Dermer (1986). 
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4.1 COMPTON SCATTERING OF THE PAIRS 

The effect of Compton scattering on tlie pair distribution can be analyzed in a manner 

similar to that of the photon Comptonization discussed in the previous section. Let p and 

Pi be the momenta of the incident electron and the photon, respectively. Let q and gi be 

their corresponding momenta after the scattering. We require the final electron energy to 

be 7. Hence, we set q = 7(1, /3), where cf3 is the velocity of the scattered electron. We write 

p = 7i(l, f3i), pi = £1(1, ki), and gi = e(l, k). Using qf = {p + Pi — qY = 0, we obtain a 

relation between the initial energy of the photon and the final energy of the electron given 

hj El = El, where 

e, = ^^^'-\ (4.1) 

ai7i - 07 

In this equation a = 1 — /3 ■ ki, ai = 1 — /3i ■ ki, and b = 1 — (3 ■ (3i. Let ^ be the cosine of 
the angle between the vectors (3 and ki. The cosine of the angle between the vectors (3 and 
f3i is defined to be //'. The angle between the planes formed by the pairs of vectors (/3, ki) 
and {f3,l3i) is defined to be 0. Now the emission coefficient due to Compton scattering can 
be written (see the appendix for details) as 



Vh) = en, (n_ + n^) r^ J d^d^'d^j, F.(7i) F,{e,] 



aiX 



IQiT E'y pi 
where 



dEi/d'j 



1 + dE/d'y 



(4,2) 



P2 Pi \Pl P2J \Pl P2J 

while pi = aiEi'-fi and p2 = aii'-f. The integration region is given by —1 < p, p' < 1, 
< < 27r, and 7min < 7i < 7max subject to the condition that Emm < £^1 < ^max- 

Next we consider the absorption coefficient due to Compton scattering. Let p = e(l, k) 
and q = 7(1, /3) be the initial momenta of the photon and the electron, respectively. The 
photon energy in the rest frame of the incident electron is given byx = p-g = 7e(l — (3fi), 
where p is the cosine of the angle between the vectors j3 and k. As in section 3.1, it can be 
shown that 

CTl f 

Xil) = ^ J dfidEF,{E) (1 -/3/i)atotai(a;), (4.4) 

where cxtotai is given by equation (|3.4|). The integration domain is given by — 1 < /i < 1 and 



- < E < E 
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4.2 PRODUCTION AND ANNIHILATION OF THE PAIRS 

The analysis for this case is analogous to that for the pair annihilation emissivity dis- 
cussed above. Let pi = £j(l,kj) be the momenta of the photons in the C-frame, where £, 
are their energies and kj are their directional unit vectors. Let p = 7(1, /3) be the mo- 
mentum of one of the emitted particles. Here c/3 is its velocity in the C-frame and 7 
is the corresponding Lorentz factor. The momentum of the C-frame itself is denoted by 
q = (1,0). We denote the CM-frame quantities with a suffix 'cm'. Let picm = ^cm(l,kcm), 
P2cm = £011(1, -kcm), Pcm = 7cm(l, /3cm), and gem = 7c(l, /^c) represent pi,P2,P, and q, re- 
spectively in the CM-frame. The velocity of the C-frame as measured in the CM-frame is c/3(, 
and 7c is the corresponding Lorentz factor. Various directional cosines are defined as follows: 
fj.,x,y, and z are the cosines of the angles between the pairs of vectors (ki,k2), (kcm,/3cm)5 
(kcm,/3c)5 and {(3^^, (3 J, respectively. The angle between the planes formed by the pairs of 
vectors {f3^,(3^^) and (/3c,kcm) is defined to be (p. We have 7cm = £^cm = \^[£i£2{^ — /^)/2], 

7c = (61+62)/ {2 Ecm), y = (£2-£l)/(2/5c7c£:cm), Z = (7c7cm-7)/^, whcrcaS A = /5c/?cm7c7cm, 

and X = yz + y/[{l — y'^){l — z'^)] cos0. We can now write (see the appendix for more details) 
the pair creation emissivity as 

•'(-'> = l6S?/,*^^ni^,fa)<i^.lii^g) . (4.5) 

1^ ' «=i \ / cm 

where the differential cross section is obtained by multiplying the one given by equation 
( ^I6D with P^^. The integration domain is given by emin < £^1,2 < £^max, — 1 < A* < 1, and 
< < 27T, subject to the condition — 1 < 2; < 1, which is equivalent to r_ < 7 < r_|., 
where r± = 7c7cm(l ± /^c/^cm)- 

For the absorption coefficient due to pair creation, consider an electron of momentum 
p = 7(1, /3) annihilating with a positron of momentum p' = 7'(1, f3'). Their common Lorentz 
factor in the CM-frame is given by 7cm = V[77'(l — P/3'fJ')/2], where fi is the cosine of the 
angle between the vectors /3 and f3'. Setting Si = 7, Ej = 7', rij = n±, 6ij = 0, dQj = 2'7idfi, 
Fj = Fe, and J-'ij = Prlrill')'^ in equation ( |2.17| ) we find 

X±(7) = ^ / d^'df^F^ii) ^^ atotai(7cm). (4.6) 

2 J 77' 

The integration is over the region 7min < 7' < 7max and — 1 < /U < 1 without any restriction. 

Here the limiting energies of the pairs are denoted by 7min and 7max- Finally cxtotai is the 
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total cross section for the pair annihilation, which is obtained by dividing the one given by 
equation ( |3.8| ) with 2[3l^ 

4.3 BREMSSTRAHLUNG COOLING RATE 

Since this process is much slower than all other reactions (roughly by a factor of a - 
the fine structure constant) we can treat it to be continuous in the energy and momentum 
(i.e., — ^ <^ 1) and use a continuity equation to describe it. At any time t, the density of 
electrons in the energy interval (7,7 + ^7) is given by neFe{'~f)d'y. Clearly ne-Fe(7)7(7) is the 
flux density of the electrons entering this interval and neFe{'y + d'~f)'j{'-f + ^7) is that due to 
the electrons leaving this interval (notice that 7 is negative in the case of electron cooling). 
The net contribution to the electron or positron kinetic equation is now given by 

- [ne(t)Fe(7,t)] = -_[n,(t)Fe(7,t)7] = C{j,t). (4.7) 

The right hand side of this equation is essentially ATTj3'~f'^{ri — xf) for the process. The cooling 
rate I7I can be written as the sum 



|7| 



/oo 
rf7'i^e(7')[^ee(7,7')+^ee-(7,7')]- (4-^ 



The cooling rates Eee,Eee, and Eep for e^-e^, e^-e^, and e^'^-proton processes, respectively, 
are given in the appendix. 

4.4 THE EFFECT OF COULOMB COLLISIONS 

Finally we analyze the effect of Bhabha and M0ller collisions (collectively termed as 
Coulomb collisions) on the electron spectrum (see e.g.. Baring 1987b or CB90 for a similar 
treatment and Dermer & Liang 1989 for a Fokker-Planck treatment of this problem). Here 
we ignore the diffusion term (which arises from the second order derivatives with respect to 
energy) that would arise in the Fokker-Planck expansion of the kinetic equation as well as 
the contribution from the pair-proton collisions (which is a much slower process). Consider 
an elastic scattering in which an electron with momentum p exchanges a momentum q with 
a target particle in the plasma which is either an electron or a positron. In both cases the 
collision cross section diverges for |q| — > and it falls off rapidly for larger values of |q|. 
We define 6 to be the angle by which the incident electron is scattered. Small values of |q| 
correspond to the small angle collisions (0 ^ 1). More precisely, |q| = |p| ^ when 6 is small. 
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Figure 2: Breinsstrahlung cooling time for (a) an electron and (b) a positron of energy 7 in 
a background thermal plasma (of electrons only) of density ric- In the former case we use 
the cooling rate E^e and we use E^.^ for the latter. Remaining cooling rates vanish in this 
particular example. The dimensionless cooling time is defined by t^. = |7|/(cneO"Th7)- Since 
the main time scale in the kinetics of the plasma is ~ (cneO"Th), tc ^ 1 at higher energies 
means that bremsstrahlung cooling is not very efficient at these energies. 
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It is well known that the divergence of the cross section for — * results in the domination 
of the relaxation process by the scattering events with small angular deflections. In many 
situations we can completely ignore the contribution from the collisions which are producing 
large angle deflections. Let /^/2 be the distance an electron has to travel in order that its 
mean-square deflection is = 7r/2 and suppose L^/2 is the distance it has to travel so that it 
is deflected by an angle of 7r/2 in a single scattering, with a probability close to unity. It 
can be shown that L^/2 — 16 7^/(457rnerg) and -^,^72/^71/2 — 2 In Ac, where 7 is the mean 
electron momentum in the background plasma. The latter ratio, in a nonrelativistic plasma, 
turns out to be 8 In Ac but the expression for Ac is different in that case. The Coulomb 
logarithm for a relativistic plasma can be shown to be In Ac = 37 + (31n7 — lnne)/2. In 
this equation rig refers to the number of electrons per cubic-centimeter. We consider only 
those plasmas for which In Ac > a few, which means that only small-momentum-transfer 
collisions are relevant. In this limit Bhabha and M0ller cross sections are equal. Therefore, 
we do not distinguish between electrons and positrons in the foregoing analysis. Consider 
two distributions f\ and /2 of electrons. The Boltzmann equation for /i can be written as a 
continuity equation in the momentum space as 

|a(p) = -|lS;(p), (4.9) 

where SI is the flux vector in the momentum space (see the appendix for its definition). 
Combining equation ( p.8|) with the above continuity equation we obtain 



[^(7) - x(7)/(7)]i = Cuil) + C,2{l), (4.10) 

where 

CUl)=^7r'crl lnAc/5 — I rf7'/?/5'7"Q(7,70, (4.11) 

while 

r\ O "I 1 

/i(7)^/s(7')-/.(7')^/i(7) /_^rf/ii?o(7,7',/i)- (4.12) 



Q(7,7') 



The derivation of equation (|4. 1 1| ) , along with the definition of the quantities involved, is given 



in the appendix. Here Cn comes from the collisions within the electrons of distribution 
/i and C12 comes from the collisions of electrons of distribution /i with the electrons of 
distribution /2. In each case we have to use the appropriate electron density in the Coulomb 
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Figure 3: Coulomb collision time for (a) a power law distribution with an index 5 = 2 
relaxing in a thermal background and (b) a power law distribution relaxing through self- 
interactions. In both cases we have used In Ac = 20. In general the time it takes to establish 
thermal equilibrium is many times that of the collision time. The spikes in these figures 
indicate that the emission and absorption rates balance at that energy, because of the form 
of Cis (see the text). 



19 



logarithm. Clearly, Cn vanishes when /i is an equilibrium distribution. In Figure 3 we 
give two examples of Coulomb relaxation. Recall that for electrons /(7) = nF(7)/(47r/57^). 
In the first example (Figure 3a) we consider a nonthermal population of density rii and 
spectrum Fi oc 7"^ for 7 > 1 interacting with a thermal background of density n2 ^ ni 
(/2 is the background distribution) and temperature 6. The electrons relax mainly through 
collisions with the background and the reaction rate is determined by C12 above (Cn is 
negligible). The dimensionless colhsion time is defined by t^ = cn2crTh/i(7)/|C'i2(7)|- In 
the second example (figure 36) we consider a power law distribution of density rii and an 
index 6, relaxing through self- interactions (there is no thermal background). In this case 

tc = C?2i(JTh/l(7)/|Cll(7)|- 

5 The computational method 

In the kinetic theory approach to nonequilibrium plasmas that we have presented in the 
preceding sections, the computational task is reduced to evaluating many collision integrals 
(for each energy bin, after each time step) quickly and efficiently, without compromising 
on the flexibility to handle many types of distribution functions. Now we explain our new 
algorithm for adaptive Monte Carlo integration which meets this demand. Our approach is 
similar to the PSD method discussed in the introduction, with the principal difference being 
that we are not using the conventional MC method (based on uniform sampling) to compute 
the transition probabilities (the collision integrals). There are also some similarities between 
our method and the LP method described in the introduction. Both methods use statistical 
weights within a Monte Carlo scheme. In the LP method, these weights are introduced in an 
ad hoc fashion, based on the energy carried by the LPs. In our method, we use probability 
weights (see below) to enhance the sampling rate in those regions where the contribution to 
the integral being evaluated is greater. But these weights (known as importance weights) are 
generated internally, through a minimal- variance prescription (see equations |5.12| and |5.14| ). 



Therefore, what we are using is a Monte Carlo method based on importance sampling. 

Being a kinetic-theory approach, our method bears a lot of resemblance to that of C92. 
Both methods discretize the kinetic equations by following the particle and photon statistics 
with reference to energy bins which are equally spaced logarithmically. Both methods have 
similar simplifying assumptions about the isotropy and the spatial uniformity of the photon 
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and pair distribution functions. Our method is more flexible than that of C92 because of 
the adaptive nature of the underlying Monte Carlo scheme (C92 uses numerical integration 
to compute the collision kernels). However, there are several differences between these two 
methods. It is straight forward to deal with anisotropic distributions in our method, ex- 
cept that some of the collision integrals given in this paper require some modifications to 
incorporate them. Moreover, we have not taken into account, the escape of photons and 
pairs from the system. We use the exact collision kernels throughout, as opposed to the 
"two-moment" approximations used by C92 or other approximation schemes (see CB90). In 
contrast to the method of C92, our time-evolution code is entirely dynamic i.e., we do not 
precompute and store any quantities, thereby obviating the need to approximate the final 
distribution functions after each time step. The integrands of some of the collision kernels 
have very narrow peaks ("integrable singularities") which are hard to evaluate through nu- 
merical methods used before (e.g., C92 and references therein). Such "singularities" can be 
easily integrated by our method. No special attention is required because the algorithm is 
adaptive - it automatically adjusts the sampling rate, iteratively, to a high value at such 
points. Now we describe our computational method in detail. 

5.1 DISCRETIZATION OF THE KINETIC EQUATIONS 

In order to simplify the analysis, we will consider only the case for which n^ = n_ = n^ 
and -^+(7) = -^-(7) = -^0(7)- This can be easily extended to the more general case. Let the 
net collision rate for the photons, due to Compton scattering, be given by 

A^{e, t) = Ane^ [r]^{£, t) - f^{e, t)xyie, t)]^^^^^^^ . (5.1) 

The corresponding collision rate due to pair annihilation and creation (ee <-^ '-ff) is denoted 
by B^{e, t). In an analogous way we define Ac(7, t) and Be{'j, t) for the corresponding collision 
rates for pairs. From the photon rate equation 

- K(t)F^(£, t)] = A^{e, t) + 5^(5, t) (5.2) 

we obtain 

[A^{e, t) + B,{e, t)] At - F^je, t)An,{t) 

/\t^{e,t) = — , [b.6) 

n^[t) + An^{t) 

where AF^{e,t) = F^{e,t + At) — F^{e,t), An^(t) = hj{t)At, and n^(t) is given below. 
Similarly we can obtain AFe(7,t). The time increment At for each time step is chosen in 
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such a way that 

d 



p.ie,t) 



^t < ///^(e, t) and 



^/e(7,t) 



At</e(7,t) (5.4) 



for all values of e and 7. In our computation we have used u = 0.1, which means that the 
maximum change in F^ or F^ in any energy bin, during any time step, is less than or equal 
to 10%. Now we determine n^ arising from the pair processes (there is no change in rig or 
n^ arising from Compton scattering). We have 

/•OO /"OO 

/ deA^{e,t) = and h^{t) = / deB^{e,t) (5.5) 

and two analogous equations for pairs. It can be shown that the positron annihilation and 
creation rates are given by 

nann(t) = ^^ f dfidjdj' F,{j,t)F,{j' ,t)(^^^a,e^,, (5.6) 

and 

en ft) r 
ncr{t) = — ^ — / dfidede'F^{e,t)F^{e',t){l - /i)a^^_,ee, (5.7) 

respectively. Here (Jee^'yy and a^^^ee are the total cross sections, which are given in the 
previous sections. For each positron annihilated or created, there will be a creation or 
annihilation, respectively, of two photons. Therefore we have 

n^{t) = 2 [raann(i) " ^cr(t)] and he{t) = hcrit) - riann(t)- (5.8) 



We remark that n_ = n^ = rie- We have verified equations ( p.5|) and ( ^.8|) in all our 
computations for time evolution, thereby ensuring the number conservation. In addition, 
we have verified the conservation of the total energy after each time step. Now we can 
use equation (|5.3|) iteratively, to obtain the time evolution of F^ and F^ from the initial 
data viz., Fe(7,0), F^{e,0), and the initial densities. We have discretized the energy {e 
and 7) with twenty energy bins per decade of energy and used a logarithmic interpolation 
between these points to reconstruct F^ and F^ for the subsequent time steps, which are then 
used in the collision integrals. Now the problem reduces to an efficient evaluation of these 
multidimensional collision integrals with complicated integrands. For this purpose we have 
developed a new version of an adaptive and iterative Monte Carlo method. It progressively 
adjusts itself to the nature of the integrand. We describe our algorithm below. 
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5.2 THE ADAPTIVE MONTE CARLO METHOD 

A general purpose algorithm for multidimensional integration which is widely used in 
the experimental particle physics is given by Lepage (1978). It is an iterative and adaptive 
scheme. A computer program implementing this method, known as VEGAS, can be found in 
Press et al. (1992). However, we have found that it has several shortcomings when applied 
to the type of integrals that arise in the kinetic theory. Not only is the convergence weak 
in some cases, we have found that the subroutine gave erroneous output for the high energy 
tails of the distributions. This is a significant obstacle because of the integrals over energy 
that we have to perform at the end of each time step. That integration makes the errors 
propagate to lower energies (where the results are otherwise accurate) during the succeeding 
time steps. We will briefly explain the original method by Lepage and then describe our 
modified scheme which can handle the integrals we need. Firstly, by scaling the integration 
variable, any multidimensional integral can be written in the form 



1 = j dvfir) 



(5.9) 



where r = (zi, 2:2, ..., -2„), dr = IliLi dzi, and the integration is over the n-dimensional hyper- 
cube < Zi < l,i = 1,2, ...,n. If we generate M ^ 1 random points r^ with a normalized 
probability density p(r) then the integral can be approximated by 

' fin) 



J 



M ^ pivk) 



(5.10) 



The variance is given by 



Ap\ 



M-1 
1 



dr 



p(r) 



I' 



M 



y /^(^^fe) _ j2 



The optimal choice for p(r) which minimizes the variance is derived from 
which implies that 



(5.11) 



p(r) 



l/(r)| 



(5.12) 



(5.13) 
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This is the central theme of the importance samphng technique - sample more in the regions 
where the absolute value of the function is larger. However, observe that the denominator is 
the integral itself! Thus we need an algorithm to solve it iteratively, starting with a reasonable 
guess for p. Then we calculate the integral by using equation ( |5.10| ) which then determines the 
new form for p{r), and so on. If this process converges in a manageable number of iterations, 
then we will have achieved our goal. The data storage requirements of directly implementing 
this scheme are well within the reach of many present-day computers. The method by Lepage 
consists of a restrictive assumption that the probability density is separable. For instance, 
when n = 2 and r = {x, y), the separability means that p{x, y) = Px{x)py{y) and to minimize 
the variance we need 

T—\'^'^[Px,Py]+K dxpx{x) + XyJ dypy{y)\ = 0, (5.14) 

which implies that 

(5.15) 



Px[X) 



1/2 • 



Io'dx[j,Uyi^ 

and a similar equation for Py{y)- For arbitrary dimensions, this scheme is implemented in 
the VEGAS subroutine, mentioned before. The motivation for assuming the separability, 
according to Lepage(1978), is that it limits the storage requirements. It is not a good 
assumption in general. Therefore we proceed to implement importance sampling directly. 
All essential features of the algorithm can be captured in a one dimensional example which 
we will consider first. Then we will show how it can be generalized to higher dimensions. 
Consider the integral I = Jq dxf{x). Let p{x) be the normalized probability density we 
want. Suppose N is an integer greater than unity and = xq < xi < X2 < ■■■ < xn = i, 
while Axj = Xj — Xj_i for i = 1,2, ..., N. We will use the following discrete representation of 
the probability density: 

P(^) = TTTi — ' if ^«-i <x<Xi, (5.16) 

JM AXi 

so that /^.'_ dxp{x) = 1/N for all i. Here the bin sizes Axi need not be all equal but all bins 
have the same probability weight. If the bin sizes are equal, we will get a uniform probability 
distribution leading to the crude Monte Carlo method. Now the integral is approximated 
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by X = Z^fcli fio,k)/Mp{ak), where < a^ < 1 are uniformly distributed random numbers. 
Typically M :§> N. Let 

«^ = T^Ec.(fc)l/K)l, (5.17) 

where Ci{k) = 1, if Xj_i < fl/t < Xj, and is zero otherwise. Clearly, Y^iLi^iAxi = I. 
Therefore Wi = MjAxj/X is the importance weight associated with the i*^ bin. Since different 
bins contribute different amounts to the integral, the idea now is to find a new set of bin 
spacings {xi,X2, ...,XAr_i}, so that all bins have equal importance weight wq = l/N. Let / 
be an integer (which depends on the bin location i) such that 

/ i+i 

Y^ W.m < iWo < Yl "^rn- (5.18) 

rra=l m=l 

Then the new grid position for the i^^ bin can be obtained from 

Xi,ne^ = xi^oid H [iwo- y2 '^rn] (a;/+i,oid " a^«,oid) ■ (5.19) 

^0 V m=l / 

However, in practice we must damp the convergence so that the contribution from the low- 
importance bins is not overly suppressed. As in the method by Lepage, we will damp the 
convergence by using the modified importance weights given by 

1 — Wi 



w' 



(5.20) 



_log(l/wi)_ 

which gives w'q = Y.f=i '^[/N . We now replace wq and Wi with the corresponding primed 
quantities in the above equations. The new probability density is now determined by using 



4, new — 



equation ( 5.16 ) and the process is repeated iteratively. If it converges, we will have x 
a^j.oid for all i, from which we can obtain the desired estimate for X. Now we give the extension 
of this scheme to two dimensions. We will assume that the number of bins is N for each 
dimension. A discrete representation of the probability density is given by 

Pi.^^y) = i;:7^ x^) if Xi-i<x<Xi and yj-i<y<yj. (5.21) 

Jy'^Axi Ayj 

This does not mean that the probability density is separable because Axj and Ayj are not 
independent in general. The integral is now estimated by 
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where < a^ < 1 and < 6^ < 1 are uniformly distributed random numbers. Let 

Ar2 ^ 
^ij = 17 Z! ^»i(^) \fi(^k, hk)\ , (5.23) 

^^^ k=i 

where Cij{k) = 1, if Xj_i < a/. < Xi and i/j^i < bk < Vj, and is zero otherwise. Now we 

define Ui = Y.f=i hijAyj and Vj = Y.f=i hijAxi. Clearly, X = Y.f=i UiAxi = Y^jLi VjAyj. Let 

Wxi = UiAxi/I and Wy^ = VjAi/j/I. From these importance weights for x and y grids we 

can obtain the corresponding damped weights and proceed to iterate as if these were two 

one-dimensional problems. Generalization to arbitrary dimensions is now trivial. 

In all our applications we found that the values A^ = 70 and a = 1.3 (for the damping 
index) gave stable and satisfactory results within at most ten iterations or so. In general it 
is advisable to start with a few thousand samples and after several iterations, increase M 
(and retaining the resulting grid) and further iterate, and so on. For many types of integrals, 
of at most five dimensions, we found that Mmax — 10^ samples to be adequate. For all the 
results presented in this paper, we have used the subroutine ran2 in Press et al.(1992) for 
the random number generation. We find that our method is faster than the crude Monte 
Carlo method (using uniform sampling) by a factor of ten or better, which is also the case 
with the method by Lepage (when it is apphcable). 

6 Time evolution and equilibria 

Here we give an analytical description of the equilibrium states of a pair plasma, in 
terms of the initial conditions. For two specific examples, we follow the relaxation toward 
equilibrium using our time-evolution code. These examples are meant to demonstrate that 
the whole formalism of this paper (the collision integrals and the computational method) 
actually works. We are considering a homogeneous, stationary, isotropic, and nonmagnetic 
system. There are no radiative transfer or hydrodynamic effects. On short time-scales 



t ^ tTh = (^+o"Thc)~^ the kinetics is determined by the rate-equations alone (see eq. |p.8|| ). 
We have seen that the collision integrals for these equations are nonlinear functionals of the 
distribution functions. Given the initial state of the plasma, we can solve these first order 
coupled and nonlinear integro-differential equations to determine the time evolution of the 
distributions. The system is characterized by the densities n^, n+, and Up and the spectra 
F^{6) and Fe['j), all of which depend on time. Their values at t = define the initial state 
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of the system. The total density of the particles is given hj n = n^ + 2n+ + Up and the 
total energy density (including the rest energy of the pairs) is given by -u = u^ + u^ + u+, 
where u^ = n^ J^ eF^{e)de and u± = n± J^ 'yFc{'y)d'y. The mean energy per particle is 
given by £ = u/n. We see that there will be no change in n due to Compton scattering or 
the pair annihilation and creation. It will change only as a result of bremsstrahlung (also 
double Compton scattering and the pair annihilation into three photons) which operate on 
a longer time scale t ~ txh/tt (a is the fine structure constant). However u remains constant 
throughout. Therefore we can divide the approach of the system toward equilibrium into 
two phases: (i) The faster phase in which both u and n remain constant and the system 
approaches to a state of kinetic equilibrium so that the total reaction rates for Compton 
scattering and the pair annihilation vanish (separately). This state is characterized by a 
temperature and the chemical potentials fi^ and jl± (ii) The slower phase in which u is 
constant but n changes, mainly due to bremsstrahlung (or its inverse, and other radiative 
processes) so that the system finally reaches a thermal equilibrium state characterized by 
a temperature Oq and a total density hq. In this state the chemical potentials vanish (see 
below). If ©0 < then no > n, which means that this phase is mainly the cooling of the 
plasma through bremsstrahlung and other similar processes. On the other hand, if Oq > <9 
then the plasma will heat up due to the inverse bremsstrahlung (free-free absorption) and 
other radiative processes. 

6.1 KINETIC EQUILIBRIUM: THE DENSITIES AND THE TEMPERATURE 

Consider Compton scattering of an electron of energy 7 and a photon of energy e. The 
respective energies after the scattering are taken to be 7' and e'. If the total reaction rate 
vanishes, then we have 



mf-rie) 



3 



1 + ^f.ie') 



f\i)f,{e') 



1 + |/.(^) 



(6.1) 



where we have retained the Bose-Einstein enhancement factor for the photons and Aq = 
h/mc. The factor half in this equation takes into account the polarization degeneracy of the 
photon states. Using the general form of the distribution functions 

/'y(^) = r 7 ^ T and /±(7) = T^exp — — (6.2) 

'^ ' Ag[exp(^)-l] '' \l \ Q± ) 
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and equation ( |6TT|) we find 9+ = 9^ = 9_. We denote this common temperature by 9. 



Notice that equation (3J.) does not yield any condition on the chemical potentials. Now 
requiring that the total reaction rate should vanish for the pair annihilation and creation as 
well, we find 



/+(7+)/-(7-) 



1 + yAl^i) 



1 + §Ue2) 



A (^0/7 (^2), 



(6.3) 



where 7-1- are the pair energies and £1^2 are the photon energies. Using the fact that the pairs 
and the photons have a common temperature 9, we obtain from this equation Jl_+fl^ = 2fl^. 
If there are no ions in the plasma (i.e., rip = 0) then n_ = n+ so that /i_ = /!+ = Jl^. By 
assuming that exp[(5 — /i^)/9] ^ 1 and exp[(7 — /i±)/9] ^ 1, for the relevant energies, we 
obtain the distribution functions in the kinetic equilibrium state to be 



The densities are given by 



V 9 



and /±(7) = T^exp 



Ag 



9 



(6.4) 



n^ 



and 



n± 



/9 

Aixe^ f^{e)de = 167r — 



8n 



exp ^^ 



9 



47r7-y/72 - l/±(7)rf7 = ^©^2 ( ^ ) exp 



9 



(6.5) 



(6.6) 



where Kn is the n^^ order modified Bessel function of the second kind. Using the relation 

2/i^ = jl_ + jlj^ we find 



n^ = A(^n^{np + «+), 



(6.7) 
n, we obtain the 



where ( = Q /K2{l/Q). Finally, from the equation n^ + 2hj^ + n^ 
densities n^ and n_|_ in terms of n and np. When (^ 7^ 1 we obtain a quadratic equation for 
fi+. It turns out that only one of its roots is physical (i.e., both n^ and n_|_ are non-negative). 
The physical root is given by 

1 
2 



n+ 



[n-n,){l-eV-n, 



(6.8) 



where n^ = (y/ln"^ — (1 — C^)?2p]. When C = 1 (which is true when 9 = 0.493) we get 
n^ = [n — UpY / {An) . Therefore we have the necessary densities in terms of the temperature. 
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When there are no ions (rip = 0) these solutions take a simple form given by 



n_ = n I 



02 + i^2(i/e) 



■n and n^ 



92 + ^2(1/6) 



-n. 



(6.9) 



Now we determine the temperature in terms of the initial data. We have 



M^ 



Aixe^ f^{e)de = 3<8)n^ 



and 



u± 



7 yi^ - l/±(7)c?7 



307^2(1/6) + iri(l/e). 



n±. 



(6.10) 



(6.11) 



'1 ■ i^2(l/e) 

Using the energy conservation equation -u = -u^ + -u^ + Uj^, we get the temperature as an 
implicit function of u, n, and Up. In the limit where rip = 0, we have 



u^ 



-n and u_ 



u. 



-n. 



92 + 7^2(1/6) 62 + 7^2(1/6) 

In this case, the equation for the temperature takes the form 

36=^ + 367^2(1/6) + 7^1(1/6) = £ [e^ + 7^2(1/6)" 



(6.12) 



(6.13) 



where e is the mean energy per particle (which is determined by the initial conditions). 

6.2 THERMAL EQUILIBRIUM: DENSITIES AND THE TEMPERATURE 

Here we determine the final temperature and densities resulting from the radiative 
processes in the second phase. We have /i_ + /i+ = 2yU^ = 0. Let /i+ = — //_ = /io and 
z = exp(/io/6o). Clearly 



^- - ^ ±1 „„^ ^ _ 167r6|; 



n± = — 3-6oi^2(l/6o)-2 and n^ 

Using the fact that n_ = Up + n+ we can show that 

167r 



X3 

^0 



n_ + n+ 



X'n 



6oi^2(l/6o)v^r+^, 



(6.14) 



(6.15) 



where x = AgUp/ [167r6o-ft'2(l/6o)]- In the nonrelativistic hmit (60 <^ 1) the pair density is 
given by 



n_ + n+ 



\3 
■^0 



{2neof'exp{-l/Qc 



^ 15^ 105^2 
1 + Y^° + 128®» 



Vl +X2. 



(6.16) 
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It can be shown that the pair energy density is given by 

IGtt 



u.+u+ = —^ 39^7^2(1/60) + 0o^i(l/eo) Vl+x2. (6.17) 

^0 



Finally, energy conservation imphes 



^^5Q4 



° + 167r 



39^7^2(1/60) + 00^^1(1/60) Vl + x^ = Xlu, (6.18) 



15 

where the first term on the left-hand side is the contribution from the photons. We can solve 
this equation for 9o in terms of u and x (equivalently rip). This completes the analytical 
description of the thermal equilibrium state in terms of the initial data. This treatment is 
exact and is valid for all energies (relativistic or otherwise) and densities (so long as the 
plasma is nondegenerate). 

6.3 TIME EVOLUTION OF THE SPECTRA: TWO EXAMPLES 

Now we consider the time evolution of the plasma for two specific initial conditions. In 
the first case the initial photon and the pair distributions are flat (i.e., F is constant) and 
nonzero within the energy (in MeV) interval 0.1 < emc^ < 10 and 0.1 < (7 — l)mc^ < 10. 
The initial densities are taken to be n^ = n+ + n— = 2 x 10^*^ cm~^. For this case we find 
a kinetic-equilibrium temperature 9 = 3.43 and the corresponding densities are found to 
be Uph = 1.36 X 10^*^ cm^ and n_ = n_|_ = 1.32 x 10^° cm^. Monte Carlo evolution of the 
spectra for this case are shown in figures 4 and 5. They agree well with the analytical kinetic- 
equilibrium solutions. For this case, as well as the second one, we have used tTh = (cfTxh^)^^, 
with n = 2 X 10^° cm^. In the second case we start with the same densities of the photons 
and pairs and the initial distributions are confined to the same band width as above. The 
only difference is that F^{e) oc e'"^ and -^0(7) oc 7"^, with suitable normalizations. In this 
case we obtain a kinetic-equilibrium temperature 9 = 0.663 and the corresponding densities 
are found to be rip^ = 1.73 x 10^'^ cm~^ and n_ = n+ = 1.1 x 10^° cm~^. Monte Carlo 
spectra for this case are shown in figures 6 and 7. Once again they are in a good agreement 
with the analytical solution. We have verified the number and energy conservation after 
each time-step. The final densities are found to agree with the predicted values within an 
accuracy of 10% or better (which can be improved by using more energy bins). In both cases 
the kinetic-equilibrium solution is moderately relativistic. It is clear from figures 5 and 7 that 
the low-energy part of the spectrum relaxes before the high-energy end. The cross sections 
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Figure 4: Time evolution of the photon spectrum (sohd hne) starting from a flat initial 
spectrum (dashed line). Initial pair spectrum is fiat as well. It is clear that the softer end of 
the spectrum relaxes first. The same phenomenon is observed in the pair distribution (not 
shown here). 
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Figure 5: Final Monte Carlo spectra (the solid histograms) at t = 45tTh compared with the 
analytical solution (dashed curves) for (a) the photons and (b) the pairs, starting from the 
flat initial spectra. See section 6.3 for details. 
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Figure 6: Evolution of the photon spectrum (sohd hne) starting from the power law {6 = 2) 
distributions of the photons (dashed line) and the pairs (evolution not shown here). As in 
the previous example, the relaxation is faster at lower energies. 
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Figure 7: Final Monte Carlo spectra (the solid histograms) at t = 40tTh compared with the 
analytical solution (dashed curves) for (a) the photons and (b) the pairs, starting from the 
power law initial spectra. It is evident that the high-energy tails persist for a long time, 
becoming steeper with time (analogous to the relaxation in a non-relativistic plasma), but 
the number of particles (and the energy) in these tails is less than a few percent of the total. 
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(and hence the reaction rates) decrease with the energy, thereby making the relaxation slower 
at higher energies. A part of the deviation from the analytical solutions that we see in figures 
5 and 7 (in the high-energy tails) could be an artifact of our sparse (logarithmic) binning 
at higher energies. It can be rectified by using more bins in the high-energy end (more 
computing time). For the above cases the final thermal-equilibrium temperatures turn out 
to be 60 = 4.36 x 10"^ and 0o = 2.98 x lO^^^ respectively. 

7 Conclusions 

We have developed a new computational method for solving the Boltzmann equa- 
tions of a pair plasma which is applicable for arbitrary energies (in the X-ray and 7-ray 
bands), densities, and distribution functions. We have fully analyzed all relevant micro- 
scopic processes in a pair plasma viz., Comptonizaion, the pair creation and annihilation, 
bremsstrahlung and the associated cooling, and Coulomb collisions. The spectra from the 
individual collision integrals, using our expressions and the numerical method (for Compton 
scattering, pair annihilation, and bremsstrahlung), are in a good agreement with several pre- 
vious results obtained by using different methods (e.g., S82a; CB90; and Dermer 1986). The 
analysis given in this paper can be very easily extended to an inhomogeneous and anisotropic 
plasma. It will only change some of the collision integrals and add a spatial component to 
the kinetic equations. That will result in an increase in the computational time but it will 
still be manageable by the present day work stations. Presence of the magnetic fields will 
alter the kinetics (through synchrotron emission) and it can be modeled along the same lines 
as that of C92. We have developed a modified version of the adaptive Monte Carlo method 
which is very efficient and robust. It is faster than the crude Monte Carlo method (using 
uniform sampling) by at least a factor of ten and is more flexible than the numerical inte- 
gration methods (which do not use random sampling) which are used in the past. We have 
obtained the analytical equilibrium solutions for a general set of initial conditions. Finally 
we have tested our Monte Carlo evolution scheme for two specific sets of initial conditions 
and found that the results compared favorably with the corresponding analytical solutions. 
The method is found to be very stable. In each of the examples considered, the program has 
analyzed a total of ~ 10^° collision events. This stability, accompanied by its generality and 
the inherent fiexibility, makes this technique suitable for many astrophysical applications. In 
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particular, this formalism can be applied to the expanding pair plasmas in the 7-ray-burst 
sources in their final stages of evolution (when they are only moderately optically thin), 
AGN, and the emission from hot accretion discs near black holes. 
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Appendix 

COMPTON SCATTERING RATE FOR PHOTONS 

The cross section in the C-frame is given by 

da 1 da ^. . , . ^. 

with 6{6 — e) = 6{ei —eijdsi/de. Here dVt is an infinitesimal solid angle around the direction 
k (similarly dVt' is defined with respect to k'). It is easy to see that dei/de = ^^ aai^. Now 

(da\ ^ dil(da\ 

\dn L , dn [dnL , ' ^ ' ' 

\ / C— irame \ / K— rraine 

Since e'^dn = e''^dVL' we get dVt/dVt' = (70)^. Finally, 

(^) - i <-) 

\ / R— frame ^ 

is the Klein-Nishina formula in our notation, where A = ^^ — ^ sin^ 9' + 1. This leads to 
In equation ( |2.13| ) we set /3rei = 1, -^12 = ai, rii = n^, n2 = n_ + n^, Fi = F^, F2 = F^, 



62 = 7, and 612 = 0. Clearly dQidQ2 = 2TTdfidfi' d(j), a = l — l3fi',ai = l — (3 fi", where 
/i" = /i/i' + ^[(1 — /i^)(l — /i'^)] COS0, and b = 1 — fi. These substitutions lead to equation 

(PD- 



Bremsstrahlung emissivity 
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Here we derive equation (|3.9| ) and explain the notation used in that connection. We 
are interested in the processes in which two particles of momenta pi = 7j(1,/3j), i = 1,2 
radiatively scatter on each other to produce a photon of momentum p = £(l,k). Here c/3j 
are the particle velocities in the C-frame and 7^ are the corresponding Lorentz factors, e is 
the energy of the emitted photon, and k is its directional unit vector. Let /i, /i', and /x" be the 
cosines of the angles between the pairs of vectors {fB^, f32), (/3i, k), and (/32, k), respectively. 
The angle between the planes formed by the pairs of vectors {j3i, [3^ and (/S^, k) is defined 
to be 0. We have /i" = ft ft' + ^/[{l — /U^)(l — /u'^)] cos0. In equation ( |2.13| ), because of 
the isotropy of the distribution functions, we can write dQidQ2 = ^iidfidQ, where dQ is an 
infinitesimal solid angle around k. We define {da/de)i = s^ j dQ{da/dP)i. The case i = 1 
refers to the e^-e^ process and the case i = 2 refers to the e^-e"^ bremsstrahlung. It is shown 
by Haug (1975b) that 

f da\ arle r ,_ d 



fdQ^, iie<e*, 
J oAi 



de J ■ 71 J pAj 

= 0, otherwise, (A. 5) 






where Ai = uVu^^^^, A2 = 2y/C - 1, and p = ^{u? - 2(xi + X2)], while uj = ^/[2 {( + 1)] 
C=Pi-P2 = 7172(1 - P1P2IJ'), xi=p-pi= £71(1 - Pi /i'), and X2=p-P2 = £^72(1 - P2 1^' 
Here a is the fine-structure constant. Finally 

C- 1 
^* = ^ (A. 6) 

71+72- V'(7i + 72)'- 2(C + 1)' 

and ^^^^_ 

Q = I^^KEI J Adn'\ . (A.7) 

The cross section Ci was computed by Haug (1975a, eq. [Al]) and C2 by Haug (1985a, 
eq.[Al]). This latter cross section has some minor errors and the corrections are given in 
Haug (1985b). For Ci 2 we have followed the notation of Haug except that dQ' was called 
dQp' in Ci and it was called dQq' in C2, in those papers. Going back to equation ( p.l3| ) 



we have to set 612 = 1 for the case i = 1. Hence nin2 -^ 2(^+ ~^ ""-)■ ^^ ^^^ second case 
612 = and nin2 -^ n^n_. We have Fi = F2 = F^. and e^ = 7, for i = 1, 2. With these 
substitutions the desired result follows. In the present notation J^u = y/C^ — I/71 72 and 
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dfl = dfi' d(j). The integration domain U is specified by 7min < 7i,2 < 7max, —1 < /W, /w' < 1, 
and < < 27r, subject to the condition that e*(7i, 72, /i) > £^- 

COMPTON SCATTERING RATE FOR PAIRS 

The cross section can be written as (see e.g., JR80) 

^2 



a 



- f dTf5^*\q + q^-p-pr)X, {kl 

)i J 



2e7pi 

where drj = d^qd^qi, q = 7/3, qi = ek, while X is given by equation ( [4.3| ). Here pi = 
p ■ Pi = q ■ qi and P2 = P ■ qi = Pi ■ q- In equation ( |A.8| ) we can remove three of the delta 



functions by integrating over d^q'. Using the conservation of three-momentum we obtain 
e"^ = (7i/3i + ^iki — 7/3)^ and de/d'-y = (7/? — 7i/?i/u' — eifi)/Pe, while /i and /i' are defined 
in section 4. After some straight forward manipulations we find 



2 



da Tg 



(iP 2-fepi 



dei/d'y 



X6{ei-e\), (A.9) 



1 + de/d-f 

where dP = d^q = jS'j'^ d'-f dQ, and dQ is the infinitesimal solid angle around the direction 
/3. Now, in equation ( p.l3| ) we set P^ei = 1, ^12 = ai, rii = n^, n2 = n^ + n^, Fi = F^, 
F2 = Fe, 612 = 0, and €2 = 71. Clearly a = 1 — (3 p, ai = 1 — Pifi", where /i" = fifi' + 
^[(1 — /U^)(l — yu'^)] COS0, and b = 1 — (3(3ifi'. Finally, dQidQ2 = 2tt dfidfi' d(j). These 



substitutions lead to equation (^3). 
Pair creation rate 

Let dQi be the infinitesimal solid angles around kj for i = 1,2. Infinitesimal solid angles 
around f3 and (3^^ are denoted by dQ and dQcm, respectively. In equation (|2.13|) we have 
dflidQ2 = 27idfidfi, because of the isotropy. We define da/d'j = P'j'^ J dQ{da/dP), where 
dP = P'y'^d'jdfl. We set rii = n2 = n^, S12 = 1, J-'u = 1 — /^, and Fi = F2 = F^. It can be 
shown that 



^ = h^-'^ (^).. = /^^^-^ (^)^^i^(ec.)5(7c. -.c.), 



(A. 10) 



where H is the Heaviside step function which is zero for negative arguments and is unity 
otherwise. The latter imposes the pair creation threshold. It can be easily seen that dVtcm = 
dz d(j). The delta function in the last equation ensures energy conservation. It can be written 
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in the form 5(7cm-£cm) = \dz/d^cm\ 5{z-z), where z = (7c7cm-7)^~^ and A = /3c/3cm7c7cm- 
This is the solution to the equation 7 = 7c7cm(l — /?c/?cm-2^) (i-e., P ■ q = Pcm ■ Q'cm)- Finally, 
\dz/d'y\ = A^^. After all these substitutions in equation (|2.13| ) we arrive at equation (^.5|). 

BREMSSTRAHLUNG COOLING FUNCTIONS 

Here we give the cooling functions used in equation (^78|) . The energy radiated per unit 
time in e^'^-proton collisions is given by 

/•7-1 ( da\ 

Eep{^) = cnpJ^ dee\ — \ , (A.ll) 

\ / proton 

where the protons are assumed to be at rest. Here e is the energy of the emitted photon 
and da/de is the cross section (see e.g., JR80). For Eee we start from equation (3.5) of Haug 
(1975b). After some algebra we arrive at 

Eeeiin') = "^''^^''"^ ^ r dfip^QUscPc), (A.12) 

ZHe 77' J-1 

where Ec = v^[(C + l)/2], Pc = ViiC ~ l)/2], and ( = 77'(1 — PP'fi), while /x is the cosine 
of the interaction angle. Averaging over this angle (/i-integration) gave rise to the factor of 
half above. Presence of rif, in the denominator is a consequence of our definition of Fee- The 
cooling function Qee, which is accurate to ~ 6% or better, is given by equation (3.15) of 
Haug (1975b). We reproduce it here for convenience: 



2 



2Pc 



Qee ^ 8ar^ 

e 



2' 



1-^- + ? f2 + §)ln(e, + p. 



(A.13) 



3 Ec 3 \ eI^ 
where a is the fine-structure constant. For e^-e^ process we get 

chmTI^ 'y -\- 'y f^ 

Eee{l,i) = ^: — / dfiPcQeeiScPc)- (A.M) 

zrig 77' J-i 

The cooling function Qeg is given by equations (26) and (28) of Haug (1985c). For the sake 
of convenience, we reproduce it here: 

r far^ Et=o (^^Pl if E, < 300 KeV, 

(5ee-=< „/ _A (A.15) 

[ IQar; [Ec In (e^ + Pc) - q£c + Ei=o h^c ') otherwise, 

where oq = 1.096, ai = -0.523, aa = 0.1436, 03 = 1.365, 04 = -0.532, 60 = -0.726, 
61 = 1.575, and 62 = -0.796. Here E^ = mc^Ec- 
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Landau collision integral for Coulomb collisions 

The flux vector (see Lifshitz & Pitaevskii 1981 - LP81 henceforth) is given by 

^i(p) = E/^V (/i(p)^/.(p') - /.(P')|-/i(p)) B^K (A.16) 

The superscripts z, j in this equation denote the components of three- vectors or tensors. In 
equation ( [A.16| ) the summation over j is imphcit. The components of momenta are given 
by p' = 7/5^ and p'' = i (3'\ for i = 1, 2, 3. We have t/^p' = p'^'^dVL'. Let C = 77' (1 - /?/?>), 
where /x is the cosine of the interaction angle. The tensor 5*-^ (see LP81) is given by 

^^' = ^^['^V^'1^3/2 [(C' - 1)^^'' - ^'^'i' - P"^"i' + i^'^" + P'niiC] , (A.17) 

where we have made some slight modifications to take into account the dimensions of the 
distributions and the momenta. This tensor satisfies the identity I]?=i -S*-' (/S-^ — 13'^) = 
0, for i = 1,2, 3. For isotropic distributions we have dfs/dpP = (3^ dfs/d'j. Using this fact 
and the previous identity we obtain 

SI = jdidn'(3'^"V,{^,^')J2B^^f3\ (A.18) 



i=i 



where 



^1(7,7') = E 



=1 L 



d d 

flil)-Q-jfsil') - fsil')g-flh) 



(A.19) 



We choose a coordinate frame in which P^ = /?, /3^'^ = 0, P'^ = /3'/i, /5'^ = P'y/Y^^JiF, and 
P'^ = 0. Also dQ' = 27Tdfi. With these substitutions we find 



^/i(7) = -/3^l 2nd^^diP(5'i'BV,{^,i), 



(A.20) 



where i? = 27r c r^ In Ac Bq and 



5n 



3^T^ (C^ - 1 - f^W - /5'%'V^ + WiiK) 



(A.21) 



77'(C2 _ 1)3/2 

The integral in equation ( |A.2CI| ) is the Landau collision integral for small angle deflections. 
This leads to the required result. 
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